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On  the  Conduction  of  Heat  in  a  Melting  Slab 

by 

Stephen  J.  Citron* 

Purdue  University 

Abstract:  A  new  method  for  the  solution  of  the  problem  of  heat  conduction 
in  a  melting  slab,  where  the  molten  material  is  immediately  removed  upon 
formation,  is  presented.  No  restrictions  are  placed  on  the  boundary  con¬ 
ditions  which  may  be  imposed  on  the  slab  and  the  material  properties  are 
allowed  to  be  temperature  dependent.  The  problem  of  determining  the  tem¬ 
perature  distribution  in  the  slab  and  the  amount  of  material  melted  is  re¬ 
duced  to  finding  the  solution  of  an  ordinary  differential  equation  on  the  amount 
of  material  melted.  This  reduction  from  a  partial  differential  equation  prob¬ 
lem  is  accomplished  by  determining  a  Taylor's  series  expansion  in  space  for 
the  temperature  distribution.  The  equation  so  obtained  for  the  determination 
of  the  amount  of  material  melted  is  of  a  form  readily  solved  numerically. 
Comparisons  with  known  results  for  a  slab  insulated  on  one  face  and  subjected 
to  a  constant  heat  input  on  the  other  face  are  given. 


♦Assistant  Professor  of  Aeronautical  and  Engineering  Sciences.  This  paper 
was  written  while  the  author  was  on  leave  at  the  Institute  of  Flight  Sciences, 
Columbia  University. 


INTRODUCTION 


Problems  involving  free  (moving)  boundaries  are  of  great  current 

interest  in  heat  conduction.  Howesver,  due  to  the  difficulties  caused 

by  the  non-linearity  of  the  problem  when  the  motion  of  the  boundary  is 

unknown  analytical  solutions  in  general  cannot  be  found.  Work  done  by 

Landau  ,  Dewey,  Schlesinger,  Sashkin,^^  and  Lotkin  among  others 

all  utilize  high  speed  computers  for  the  solution  of  the  partial 

differential  equations  involved.  Because  of  the  inherent  complexity 

of  obtaining  numerical  solutions  to  the  non-linear  partial  differential 

equations  it  is  of  interest  to  consider  if  the  problem  can  be  further 

reduced  before  any  numerical  work  is  begun. 

f  1+1 

Recently  Boley  developed  a  method  in  which  the  problem  is 

reformulated  so  as  to  require  the  solution  of  two  ordinary  integro- 

differential  equations.  By  expanding  the  equations  in  powers  of  the 

time  after  melting  starts,  an  exact  analytical  solution  is  obtained 

rs] 

which  is  particularly  useful  for  small  times.  The  author  devised 

a  method  of  successive  approximations  by  means  of  which  a  solution  may  be 
obtained  by  solving  an  ordinary  differential  equation  on  the  amount  of 
material  melted. 

The  purpose  of  the  present  work  is  to  present  a  new  method  of 

reducing  the  problem  of  determining  the  temperature  distribution  and 

the  amount  of  material  melted  to  one  requiring  the  solution  of  an  ordinar 

fql 

differential  equation.  It  differs  from  the  method  of  reference  in 
that  it  eliminates  the  need  for  the  successive  approximation  procedure 
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used  there.  The  main  feature  of  the  present  method  is  the  expansion 
in  space  of  the  temperature  distribution  in  a  Taylor  series.  The  one 
boundary  condition  on  the  temperature  which  is  not  identically  satisfied 
by  the  expansion  yields  an  ordinary  differential  equation  for  the 
determination  of  the  amount  of  material  melted.  The  initial  conditions 
for  this  new  equation  are  found  by  using  the  initial  temperature  dis¬ 
tribution  or  an  analytical  starting  solution  as  will  be  discussed  later. 

Once  this  equation  has  been  solved  (and  hence  the  amount  of  material 
melted  as  a  function  of  time  is  known)  the  temperature  distribution  at 
any  time  can  be  immediately  calculated. 

A  comparison  of  solutions  obtained  by  this  method  with  certain 
known  results  for  a  slab  insulated  on  one  face  and  subjected  to  a 
constant  heat  input  on  the  other  are  given.  The  example  used  to  demonstrate 
the  method  was  chosen  as  it  is  an  important  limiting  case  encountered 
in  ablation  problems.  In  addition,  this  problem  has  been  studied  by 
others  using  techniques  more  complex  than  the  method  to  be  presented. 

The  problem  thus  furnishes,  by  comparison,  an  illustration  of  the 
simplicity  of  the  method  developed  here. 

PROBLEM  FORMULATION 

Consider  the  slab  shown  in  Fig.  (l )  with  temperature-dependent 
thermal  properties  subjected  to  an  arbitrary  heat  input  Q(t)  on  one  face. 

No  restrictions  are  placed  on  the  boundary  conditions  which  may  be 
prescribed  on  the  other  face  of  the  slab.  Once  melting  occurs  at  the 
heated  side  the  problem  to  be  solved  requires  the  determination  of  the 
temperature  distribution  T(x,t)  in  the  slab  and  the  amount  of  material 
melted  as  a  function  of  time  s(t).  The  molten  material  is  taken  to  be 
immediately  removed  upon  formation. 
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'f  =  L  *  * 


(3) 


The  slab  is  now  of  fixed  unit  length  while T  measures  the  non-dimensional 
time  after  melting  has  begun.  The  following  dimensionless  quantities  are 
also  introduced* 


T  ¥ 

l 


Qm-  Qjrt} 
3* 

-  £  /SlrN 


*<T*>  (4) 

C  CT>  -  CCT^ 

car*'. 


S(y>  -  Ci-aS/(t3 

Eq.  (l )  then  becomes 


H  -  3l‘4  v  -  filin’* 

L  ^  6><» 


D  =c$2  Vc&aSS  o<^<t 

V  SK  ^ 


under  the  conditions  corresponding  to  those  of  Ea-  (2)  given  below 


a')  -  “To(V 

0(\,y)  -  \  (6) 

c")  3f>0,r>  --  SQ  -  V'* 

^  v  2.  ^ 

<n  G(exrt>&  .. ^  *  vr>1  **  }  mo 

’*v  * 

e")  S  (o^  =  o 

It  may  also  be  noted  from  Eq.  (uc )  that  continuity  of 
at  T  =  o  requires  that  S  =  o  for  M  ^  °° 


*  The  quantity  Qq  is  a  reference  value  of  the  heat  imput. 
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PROBLEM  SOLUTION 

"Hie  method  of  solution  depends  on  the  assumption  that  it  is  possible 
to  express  the  temperature  distribution  in  the  slab  in  a  power  series 
expansion  in  space  about  the  melting  face  z  *  1.  This  assumption  is  to 
be  checked  a  posteriori  by  verifying  that  all  equations  and  conditions  are 
satisfied  by  a  convergent  series  or,  lacking  that,  that  the  expansion 
formally  satisfies  all  differential  equations  and  conditions  and  in  addition 
checks  some  known  results. 

The  Taylor  series  expansion  of  the  temperature  is  given  formally  by 

-  <90, ' r}  *  2j9 (i yHV '•>  *  (wMv o'-*.. •  (7) 

*1  %  hi 

where  the  coefficients  of  the  powers  of  (z  -  1 )  must  be  determined.  These 
coefficients  are  found  by  utilizing  the  boundary  conditions  at  z  «=  1  and 
the  differential  equation  governing  the  problem  Eq.  (5).  From  F.q.  (6  b,c) 
we  have 

6o,*o  -  t 

%  v  2  h 

The  higher  order  coefficients  are  then  obtained  by  using  these  expressions  with 
Eq.  (5).  For  this  purpose  we  rewrite  Eq.  (5)  in  the  form 

C  £.^0(9  *  c  (S.S'T.-O^h  ij9 

I  ^  =>V 


(9) 
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Thus  the  second  derivative  of  £  with  respect  to  z, 
evaluated  at  z  =  1,  is  given  in  terms  of  the  first  derivative  of 
0  with  respect  to  z  and  of  the  value  of  0  at  z  =  1,  both  known 
quantities  from  Eq.  (8).  In  a  similar  manner  all  higher  derivatives 
of  Q  may  be  found  from  Eq.  (9)  by  successively  differentiating  both 
sides.  The  n'th  derivative  of  &  with  respect  to  z  will  be  a  function 
of  the  n  -  1  preceeding  derivatives  which  will  have  already  been 
evaluated.  As  examples  the  third  and  fourth  differential  coefficients 
are  given  below.* 


The  term  ^ 


^  occurring  in  the  expression  for 

above  is  found  by  differentiating  Eq.  (6c)  with 


respect  to  t. 


*  It  has  been  assumed  that  it  is  possible  to  interchange  the  order  of  t 
and  z  differentiations.  Note  also  that^(l,x)  =  1  from  its  definition. 


-7- 


Once  the  coefficient*  of  the  expansion  have  been  obtained, 
substitution  of  Eq.  (7)  into  Eq.  (6d )  [the  only  boundary  condition 
which  has  not  been  identically  satisfied  by  the  expansion]  yields 
an  ordinary  differential  equation  for  the  determination  of  &  (x). 

If,  for  example,  the  slab  were  insulated  at  z  »  o  then  Eq.  (6d) 
becomes  (o,x)  *  o  and  substitution  of  Eq.  (7)  into  this  relation 

yields  the  equation 

o  -  U/r'k  -  x  •••  /n\ 

for  the  determination  of  (x). 

Note  that  the  first  term  is  the  expansion  of  0  (z,x)  is  a 
constant,  the  second  and  third  terms  contain  jS*  (x)  and  !>  (x),  the 
fourth  and  fifth  terms  contain  (x),  Jo  (x)  and  (x).  In  general 
the  In  and  +  1  terms  of  the  expansion  will  contain  S(x)  and  its 
first  n  derivatives.  It  may  also  be  seen  from  the  expansion  that  the 
highest  derivative  of  *  (r)  occurring  always  appears  in  a  linear  manner. 

INITIAL  CONDITIONS 

In  practice  the  series  for  the  temperature  will  be  terminated 
after  a  finite  number  of  terms,  say  2n  or  in  +  1.  The  n'th  order 
differential  equation  on  (x)  resulting  from  satisfying  Eq.  (6d)  then 
requires  n  initial  conditions  of  which  only  two  are  already  known, 
namely  &  (o )  =  £  (o)  *  o.  Two  ways  are  now  suggested  of  determining 
the  additional  initial  conditions  required  to  start  the  solution.  The 
two  methods  given  here  will  later  be  compared  through  examples. 
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1)  If  an  analytical  solution  for  £  (t)  is  known  in  the 
neighborhood  of  t  =  o  then  (t  )  and  the  first  n  -  1  derivatives 
of  j£>  ( t)  may  be  calculated  at  some  time  t|  within  the  range  of 


validity  of  the  solution.  The  values  of  £ CO ,  c re*  ...  cO 

thus  calculated  are  used  as  initial  conditions  for  starting  the 

numerical  solution  of  Eq.  (6d )  at  t  =  •  The  solution  proceeds 

by  substituting  these  values  into  Eq.  (6d )  from  which  a  value  of 

^ is  calculated.  Having  a  value  of  ^ 

*5>*r,*v  ^ 

the  values  of  jgf  (t2),  \_£  (r2),  .  .  .  ^  n  "  1  f>  (t2)  where 

'*rr 

T0  «=  x.  +  £  t  are  determined  from 


J&LTx')  -  tfi'tx')  *  W 
^  3rr  V*v 


(12) 


-^T'a~’  o>TA'"*  VY"*' 


The  process  is  then  repeated  so  as  to  cover  the  time  interval  of  interest. 
It  should  be  noted  that  the  initial  condition  on  the  temperature  Eq.  (6a) 
is  automatically  satisfied  by  using  the  exact  starting  solution. 

A  method  of  obtaining  the  analytical  starting  solution  required 

[4] 

by  the  above  procedure  has  been  given  by  Boley  who  explicitly  gives 


as  an  example  the  starting  solution  for  the  case  of  a  semi-infinite  solid 
of  constant  thermal  properties  under  constant  heat  flux  at  the  melting 
boundary.  The  analogous  formula  for  the  slab  obtained  by  the  method  of 
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reference  [4]  is  later  given  in  the  present  paper  for  use  as  a  starting 
solution  in  the  examples  which  follow.  Boley's  method,  presented  here 
in  brief  for  completeness,  consists  in  dealing  mathematically  with  a 
fictitious  solid  of  constant  dimensions  under  an  equivalent  unknown 
heat  input.  This  results  in  two  ordinary  integro-dif f erential  equations 
for  (t  )  and  the  fictitious  heat  input  which  must  be  solved  simul¬ 
taneously.  By  expanding  these  quantities  about  r  *  o  in  the  proper 
manner  it  is  possible  to  determine  an  analytical  solution  valid  about 
x  =  o. 

2)  When  the  thermal  properties  of  the  slab  are  variable  it 
may  prove  difficult  to  determine  an  analytical  starting  solution, for 
the  period  after  melting  has  begun, as  described  above.  We  therefore 
require  a  different  method  of  determining  values  of  >  ***  ^ — - 

at  t  *  o  to  start  the  numerical  solution  of  Eq.  (6d )  and  we  do  this 
by  satisfying  the  one  remaining  condition,  the  initial  condition  on 
the  temperature  Eq.  (6a),  in  an  approximate  manner.  If,  for  example, 

2n  or  2n  +  1  terms  are  retained  in  the  expansion  of  6  given  by  Eq.  (7) 
then  there  are  n  initial  conditions  to  be  determined,  namely  the  values 
of  %  j  .  .  •  •  c>  .  we  have  already  shown, 

■»r  ■£*r‘*'* 

.  * 

however,  that  jS  (0)  =  £  (0 )  ■«  o  which  leaves  n  -  2  initial  values  to 
be  found  «>*  ^  } . V*  S 

Sr* 


These  n  -  2 


?  e/v 
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values  are  used  to  match  the  initial  temperature  distribution  at  n  -  2 
points  therefore  approximately  satisfying  the  only  remaining  condition.* 
The  numerical  integration  can  now  proceed  exactly  as  given  under  method 
1  above. 

EXAMPLE 

The  general  method  of  solution  given  above  and  the  two  methods 
of  starting  the  solution  will  now  be  illustrated  by  considering  the 
problem  of  a  slab  of  constant  thermal  properties  insulated  at  x  =  1 
and  subjected  to  a  constant  heat  flux  Q  on  the  other  side.  The  accuracy 
of  the  results  will  be  demonstrated  by  comparison  with  the  conditions 
which  are  known  exactly  for  the  problem  .  These  conditions 

are; 

1)  the  total  time  for  complete  melting  given  by 

±  i 

^  2  iA  ‘ 

2)  the  melting  rate  when  s  =  9.\  £  *  1 

k.  =  m 

4  (lit) 


3)  the  temperature  distribution  at  the  start  of  melting 
and  the  temperature  distribution  in  the  neighborhood  of  this  time  derived 
by  the  method  of  reference  [4], 

*  Note  that  since  Eq.  (6b,  c,d)  will  be  satisfied  the  approximate  temper¬ 
ature  distribution  at  t  ■  o  actually  matches  the  exact  distribution  at 
n  -  1  points  and  has  in  addition  the  correct  slope  at  z  =  1  and  the 
correct  condition  at  z  >  o. 
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given  by 


The  temperature  distribution  at  the  start  of  melting  is 

[6] 


12/H4  i  )* 

+  Urn-t  )<A~^  ^ 

where  ^ 

. A 

Since  0  (o,  )  -  1  we  find  for  the  relation  between  V  0/»A  ot 


V  -  _ L 


yz  «c 


(  \  -v  2  fr r  / 2mcO 


In  the  numerical  work  which  follows  we  shall  use  the 


numerical  values 


n  *  2 


v  *  2 


Substitution  of  this  value  of  V  into  Eq.  (l6)  then  yields  *.39 
We  proceed  by  first  giving  the  analytical  solution  for 
£  (t)  valid  in  the  neighborhood  of  t  =  o  which  is  derived  by  the 
method  of  reference  [4] 


-  fco'YVa’  •*  *  V>a.  'Y  + 


(18) 
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where 


bo  “  -J1  £  do 

ifir  >r 

b  v  *  -  J-  £  On 

V- 


.HmV 


W  2.  “  2 fr'N  j  TTCLo  Wo  —  2  )fi r  CL2.  +  ^  1 

5ttv  l  '  ) 

and 


aQ  -  -2o( 


^2  +  £  £  >»<*  -v  ii  (19> 


<X»  “  1  UrT  \jo 

H 

d° 

^•2  *  <  n/tF AoVo  2ol  -v  'Hr  "t  H«A  fiT  ^  £/y\<K  &^c2  /n«t 

C  fi.  __  u  Jk-j  x 


,  V  7  *  'j 

-  e  1-  may 


The  temperature  distribution  derived  by  the  same  method 
of  reference  [W]  is  given  by 


(-T  2!  Vt  1Z1  t-Vr-MU  54. 'I 
T*  vff  C  2  •*  ziTr  2  "  2|T?  } 


^  * 

♦a.  f  2  If?  (T \  <W  )  a^f.  */&.  -  2 e 
V*  (  6  *i*  l 


W2[3t  *  U1 

v-^  1  « 


+  Xl_  1  trJU.  *_/£r 

32  iV  ®2iFf 


•+ 
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-KX  -v 

“  'T^lJrT  t  §  */jBL  ♦  J-  1  £  4*C'T*  t 

*  \6  t-Vz  J 

+  VHT4  VU*  |  ^8^  */jl  (20) 

V  i  [MT+  '4*3‘4 

«* 

+  /  I  4A>i^x  <  2m  4  */l  >  V  AinJfyt.  I2n~  */£\  ~0 

'  [Mr  +  ‘/^T4  +  !£i  y*  ) 


The  expansion  for  the  temperature  distribution  in  this  example 
obtained  by  retaining  the  first  five  terms  of  Eq.  (7)  is  given  by 


8(^/0-  \  *  i<\-!LV£:c  SMyx'i  li-lVifUvO1 

v  a  n  d  *7  in  ^-j- 

4  (  <»  -  TT.^  ViQ 

L  y  an  2rA  3  (2i) 

4^SHU-E4r  in  t  $MSWV^](ytf 

(v-an  n  a.n^ 

Utilizing  this  expansion  the  condition  that  the  slab  be 


insulated  at  x  =  1  (z  =  o)  becomes 


-siW-  ?£ 


hi 


TT  *■  V 


6  (•  sn  -Ci 

V'>v  } 
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The  analytical  starting  solution  of  Eq.  (18)  was  used  to 
determine  values  of  jf(t)  and  $  (t)  at  t  ■  .04,  the  point  chosen  to 
begin  the  numerical  integration.  The  value  of  O  (.04)  was  then  calculated  from 
Eq.  (22)  when  first  four  and  then  five  terms  of  the  temperature  expansion 
were  retained.  Substitution  of  these  values  Into  Eq.  (22)  retaining  the 
corresponding  number  of  terms  then  yields  expressions  for  the  temperature 
distribution  at  r  <■  .04.  The  temperature  distributions  at  r  ■  .o4 
calculated  from  Eq.  (20)  by  retaining  first  four  and  then  five  terms  In 
the  expansion  are  compared  in  Fig.  (2)  with  the  exact  temperature  distri¬ 
bution  at  this  time  given  by  Eq.  (20 ).  The  figure  shows  the  good  agreement 
obtained. 

The  solutions  for  jS  (t  )  calculated  by  numerically  Integrating 
Eq.  (22)  after  r  »  .04  are  shown  in  Fig.  (3).  The  error  In  the  value 
of  Is  only  4.9$  In  the  four-term  expansion  and  Is  reduced  to  2.5$  In 
the  five-term  expansion.  It  may  also  be  seen  from  Eq.  (22)  that  the 
correct  final  melting  rate  »  2  FA  Is  reached  In  both  the  four 

^  IL 

TT  V 

and  five-term  expansions. 

It  remains  to  illustrate  the  method  of  solution  when  an 
analytical  starting  solution  Is  not  known.  For  this  purpose  seven  terms 
were  retained  In  the  expansion  for  the  temperature.  The  expansion  now 
contains  terms  up  to  S  (t )  which  enables  the  matching  of  one  additional 
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point  of  the  approximate  temperature  distribution  with  the  exact 
distribution  at  t  =  o  as  was  previously  discussed.  The  point  chosen 
was  the  insulated  fact  X  *  (r  *  o).  A  comparison  of  the  approximate 
distribution  thus  calculated  with  the  exact  distribution  given  by 
Eq.  (15)  is  shown  in  Fig.  (!*■)•  The  agreement  between  the  two  temperature 
distributions  is  again  good. 

The  results  of  numerically  integrating  Eq.  (22)  starting  at 
-r  =  o  with  the  additional  terms  retained  is  shown  in  Fig.  (5).  A 
comparison  is  made  in  the  figure  with  the  solution  obtained  by  retaining 
only  five  terms  of  the  expansion  and  starting  the  numerical  integration 
at  t  •  .OU  with  the  exact  starting  solution. 
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